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TECHNICAL  SUMMARY 

A  non-linear  method  is  developed  for  the  simultaneous  inversion  of 
Lg  source  spectra  and  path  Q.  The  method  is  applicable  to  both  earth¬ 
quakes  and  explosions.  The  unknowns  are  the  source  corner  frequency 
and  seismic  moment,  as  well  as  Lg  Q0  and  tj  values  (Lg  Q  at  1  Hz  and  its 
power-law  frequency  dependence)  along  paths  to  the  recording  stations. 
The  inverse  algorithm  determines  different  Q0  and  r\  values  for  each  of 
the  several  paths,  it  is  event-based,  and  it  is  suitable  for  real-time  process¬ 
ing.  The  algorithm  is  further  characterized  by  an  exhaustive,  non-linear 
search  for  source  parameters  and  linear  searches  for  Q0  and  t]  values.  The 
search  process  is  implemented  in  successive,  increasingly  refined  steps, 
designed  to  accelerate  the  computation.  With  these  characteristics,  the 
inverse  algorithm  requires  neither  a  starting  model,  nor  a  priori  informa¬ 
tion  on  the  source  and  Q;  it  thus  avoids  convergence  to  local  minima,  and 
is  fast  enough  to  be  mini-computer  based. 

The  method  is  applied  to  three  underground  nuclear  explosions  in 
eastern  Kazakhstan.  The  resulting  seismic  moment  (M0)  and  corner  fre¬ 
quency  (fc)  for  the  largest  (JVE)  event  are  1.3  (±  0.1)  x  1023  dyne-cm  and 
0.56  (±  0.02)  Hz,  respectively.  These  values  are  consistent  with  previous 
estimates  using  other  methods  and  phases,  as  well  as  the  source  scaling 
relationships  developed  for  NTS  explosions.  Based  on  the  same  scaling 
relationships,  the  M0  values  obtained  for  the  smaller  events  in  this  study 
suggest  seismic  yields  that  are  consistent,  within  reasonable  uncertainties, 
with  those  predicted  by  using  the  ISC  mb  values  with  an  mb-yield  scaling 
relationship.  The  path-averaged  apparent  Q0  and  ti  values  agree  with 
previous  estimates  in  roughly  the  same  area,  but  the  values  from  indivi- 
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dual  paths  can  vary  significantly.  The  apparent  Q0  values  are  higher  than 
750  in  two  northeasterly  directions,  and  are  between  475  and  591  to  three 
other  directions.  These  values  correlate  with  the  regional  tectonic 
features.  These  values  also  agree,  within  the  uncertainties,  with  estimates 
based  on  a  tomographic  Lg  coda  Q0  map  by  Xie  and  Mitchell  (1991). 
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Simultaneous  Inversion  for  Source  Spectrum  and  Path  Q  Using  Lg 
with  Application  to  Three  Semipalatinsk  Explosions 


ABSTRACT 

A  non-linear  method  is  developed  for  the  simultaneous  inversion  of  Lg 
source  spectra  and  path  Q.  The  method  is  applicable  to  both  earthquakes  and 
explosions.  The  unknowns  are  the  source  comer  frequency  and  seismic 
moment,  as  well  as  Lg  Q0  and  T|  values  (Lg  Q  at  1  Hz  and  its  power-law  fre¬ 
quency  dependence)  along  paths  to  the  recording  stations.  The  inverse  algo¬ 
rithm  determines  different  Q0  and  ti  values  for  each  of  the  several  paths,  it  is 
event-based,  and  it  is  suitable  for  real-time  processing.  The  algorithm  is  further 
characterized  by  an  exhaustive,  non-linear  search  for  source  parameters  and 
linear  searches  for  Q0  and  t)  values.  The  search  process  is  implemented  in  suc¬ 
cessive,  increasingly  refined  steps,  designed  to  accelerate  the  computation.  With 
these  characteristics,  the  inverse  algorithm  requires  neither  a  starting  model,  nor 
a  priori  information  on  the  source  and  Q;  it  thus  avoids  convergence  to  local 
minima,  and  is  fast  enough  to  be  mini-computer  based. 

The  method  is  applied  to  three  underground  nuclear  explosions  in  eastern 
Kazakhstan.  The  resulting  seismic  moment  (M0)  and  comer  frequency  (/f)  for 
the  largest  (JVE)  event  are  1.3  (±  0.1)  x  10a  dyne-cm  and  0.56  (±  0.02)  Hz, 
respectively.  These  values  are  consistent  with  previous  estimates  using  other 
methods  and  phases,  as  well  as  the  source  scaling  relationships  developed  for 
NTS  explosions.  Based  on  the  same  scaling  relationships,  the  M0  values 
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obtained  for  the  smaller  events  in  this  study  suggest  seismic  yields  that  are  con¬ 
sistent,  within  reasonable  uncertainties,  with  those  predicted  by  using  the  ISC 
values  and  an  mb  -yield  scaling  relationship.  The  path-averaged  apparent  Q0 
and  t)  values  agree  with  previous  estimates  in  roughly  the  same  area,  but  the 
values  from  individual  paths  can  vary  significantly.  The  apparent  Q0  values  are 
higher  than  750  in  two  northeasterly  directions,  and  are  between  475  and  591  to 
three  other  directions.  These  values  correlate  with  the  regional  tectonic  features. 
These  values  also  agree,  within  the  uncertainties,  with  estimates  based  on  a 
tomographic  Lg  coda  Q0  map  by  Xie  and  Mitchell  (1991). 

Abbreviated  title:  Simultaneous  Lg  Spectral  Inversion 
Keywords:  Non-linear  optimization,  Lg,  Source  Spectra,  Q,  Nuclear  Explosions 
1  INTRODUCTION 

The  seismic  Lg  wave  train  is  typically  the  most  prominent  phase  on  short-period 
records  observed  in  stable  continental  areas.  It  can  be  successfully  modeled  as  a  sum  of 
higher-mode  surface  waves  (Knopoff  et  al.,  1973),  or  a  sum  of  super-critically  reflected 
waves  in  the  crust  (Bouchon,  1982).  Because  Lg  is  so  often  such  a  prominent  phase  and 
because  of  its  stable  amplitudes,  it  has  been  used  to  measure  sizes  of  seismic  sources  via  a 
magnitude  scale,  (Nuttli,  1973).  The  methodology  employed  in  developing  the  mhLi 
scale  consisted  of  simultaneous  determination'  of  source  sizes  and  path  attenuation  using 
time-domain  data.  This  methodology  was  later  extended  to  digital  data  in  the  spectral 
domain  (Street  et  al.,  1975;  Campillo  et  al.,  1985;  Sereno  et  al.,  1988).  The  outcome  of  these 
studies,  in  terms  of  source  characteristics  and  path  attenuation,  have  provided  some  of  the 
most  robust  estimates  of  sizes  of  the  intra-plate  earthquakes  (Nuttli,  1973,  Chow  et  al., 
1980,  Alexander,  1985)  and  of  underground  nuclear  explosions  (Nuttli,  1986a,  1986b,  1988; 
Henson  et  al.,  1990;  Jih,  1992).  Lg  attenuation  has  also  been  found  to  correlate  with  the 
tectonic/structural  environments  (eg.,  Mitchell,  1981;  Kennett  et  al.,  1985;  Campillo,  1987). 
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Most  previous  simultaneous  spectral  inversions  of  source-to-Lg  radiation  and  path 
have  been  conducted  under  the  assumption  of  a  constant  Q  model  which  is  invariant  over  a 
large  area  covered  by  paths  to  all  stations  used  in  the  study.  In  addition  they  were  based 
on  large  data  sets,  typically  narrow-band,  collected  from  dense  networks  (eg.,  Street  et  al., 
1975;  Sereno  et  al.,  1988).  Campillo  et  al.  (1985)  and  Campillo  (1987)  separately  measured 
the  path-averaged  Qif ,  source-to-Lg  radiations,  and  path-variable  QLf  in  central  France 
using  a  dense  seismic  network.  The  result  indicates  that  is  highly  variable  within  an 
area  whose  dimension  is  less  than  a  few  hundred  km.  Time  domain  measurement  of 
of  underground  nuclear  explosions  (Nuttli,  1986a,b,  1988;  Jih,  1992)  also  required  significant 
variations  in  QLg  among  paths  with  lengths  of  a  few  hundred  to  a  few  thousand  km. 

With  the  deployment  of  modem  broad-band,  high  dynamic  range,  digital  seismic  sta¬ 
tions,  such  as  the  GDSN,  IRIS,  GEOSCOPE  and  CDSN  stations  on  major  continents,  high 
quality  data  has  become  available  for  real  time,  event-based  processing.  This  data  should 
permit  simultaneous  determinations  of  source-to-Lg  radiation  and  path-variable  QLg  values 
without  a  priori  knowledge  on  either  source  or  Q  values.  In  this  paper,  I  will  present  a 
non-linear  spectral  method  for  such  an  inversion,  and  will  apply  it  to  Lg  from  three  under¬ 
ground  nuclear  explosions  in  eastern  Kazakhstan  recorded  by  broad-band  IRIS  and  CDSN 
stations.  These  applications  simultaneously  provide,  for  the  first  time,  source  spectral 
parameters  and  path-variable  Q  parameters.  These  will  be  compared  with  a  prion 
knowledge,  including  available  previous  estimates  of  the  yield  of  the  explosions  using  dif¬ 
ferent  methods  and  phases,  various  source  scaling  relationships  for  underground  explo¬ 
sions,  and  a  tomographic  map  of  Lg  coda  Q  (Xie  and  Mitchell,  1991).  It  will  be  shown  that 
the  results  of  the  applications,  while  requiring  no  a  priori  knowledge  of  source  parameters 
or  Q,  are  generally  highly  consistent  with  available  a  priori  knowledge,  thus  demonstrating 
the  reliability  and  the  power  of  the  inverse  algorithm. 


2.  STOCHASTIC  MODELING  OF  Lg 


-4- 


Street  et  al.  (1975)  were  the  first  to  empirically  relate  the  Lg  ground  motion  spectra  to 
source  spectra.  Subsequent  developments  have  incorporated  the  effects  of  finite  Q 
(Herrmann  and  Kijko,  1983;  Campillo  et  al.,  1985),  and  the  stochastic  nature  of  Lg  (Xie  and 
Nuttli,  1988;  Ou  and  Herrmann,  1990;  Xie  and  Mitchell,  1990)  into  the  formulation.  The  for¬ 
mulation  adopted  for  the  inverse  method  of  this  study  is  that  given  by  Xie  and  Nuttli  (1988) 
and  Xie  and  Mitchell  (1990)  and  is  expressed  as 


A?(J)*S(f,  6i)G(Ai)exp 


«/T  i 

Qi(f) 


(i) 


where  /  is  frequency,  Atig(f )  is  the  Lg  displacement  spectrum  at  the  ith  station,  8,  ,  A,,  t, 
and  Qj(f)  are  the  azimuth,  distance,  mean  Lg  travel  time  and  Lg  quality  factor  from  the 
source  to  the  ith  station.  S (/  ,8i )  is  the  source-to-Lg  spectral  radiation  in  the  direction  of  8;, 
or  "Lg  source  spectrum".  G(A()  is  the  geometrical  spreading  factor  and  typically  has  the 
form 


C( Ai)  =  (M,rW  (2) 

where  A0  is  a  reference  distance.  Street  et  al.  (1975)  and  Herrmann  and  Kijko  (1983)  used  a 
A0  of  100  km  and  the  same  will  be  used  in  this  study.  The  reason  for  this  value  is  discussed 
in  a  later  section.  X,(f)  and  rt{f)  in  equation  (1)  are,  respectively,  the  site  response  and  the 
randomness  due  to  propagation.  In  a  laterally  homogeneous  medium,  r,  mainly  represents 
the  random  effects  of  multiple  ray  (or  mode)  arrivals  comprising  Lg,  whereas  in  a  laterally 
inhomogeneous  medium  it  also  represents  the  effects  of  3D  heterogeneity. 

For  earthquake  sources,  S  (f  ,9,  )  can  be  represented  by  the  omega-square  model  (Street 
etal.,  1975)- 


M0R(ei) 


4irpo,3  1 +f2/f? 


(3) 


where  M0  and  fc  are  the  seismic  moment  and  comer  frequency  of  the  source,  p  and  vt  are 
average  density  and  shear  wave  velocity  in  the  crust,  with  typical  values  of  2.7  gtm3  and 
3.5  km/s,  respectively.  R(0,)  is  a  combination  of  radiation  pattern  and  focusing/defocusing 
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effects  in  3D  media.  For  explosive  sources,  Sereno  et  al.  (1988)  proposed  a  simplified 
Mueller-Murphy  model: 

MoR(0I.)  i 

'•i)  -  — r  i ^  '  <4> 

4lTpo«  [l+(l-2  V)f2/fe2+f?f*/fc*] 

where  we  have  assumed  the  explosion  efficiency  coefficient  (k  in  equation  (7)  of  Sereno  et 
al.,  1988)  to  be  unity,  quantifies  the  overshoot  effect  and  is  defined  by  Mueller  and  Mur¬ 
phy  (1971): 


P 


I 

source 


where  "  I  source  "  denotes  quantities  evaluated  in  the  source  zone. 


(5) 


Simultaneously  inverting  for  fe,  M0/  Q, (/ ),  R (6,)  and  X,(/ )  is  normally  not  practical 
due  to  the  trade-offs  of  Q,(/),  fi(6,)  and  X;(/).  We  therefore  make  a  simplification  by 
assuming 


J?(0,)exp 


r 

rr/T; 

rr/T( 

- 

X;(/)  =»  exp 

- 

Qi(/)  J 

Q',(/)J 

(6) 


where  Q',  (/ )  is  apparent  Lg  Q  and  is  assumed  to  have  a  power-law  frequency  dependence 


Q'i(f)*Q'oif  ‘  (7) 

where  Q  and  V»  are  apparent  Lg  Q  to  the  ith  station  at  1Hz  and  its  power-law  frequency 
dependence,  respectively. 

Substituting  equations  (2),  (6)  and  (7)  into  (1),  we  obtain 


A  if)  »  ' _ exp 

VdA 


,W.- 

*1 _ 

Q'otV)  i 


ri(f) 


where  S  (f )  is  given  by 


S'if) 


M0  i 
lirpv,3  1  +/2/// 


(8) 


(9) 
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for  earthquakes,  arid 

_  M0  j 

S^(/ }  -  — •  - -  (10) 

4*po»  [i+(1-2P)/2//c2+pV4//c4] 

for  explosions.  Q  and  rj'j  in  equations  (6)  through  (8)  will  be  affected  not  only  by  the  Lg 
quality  factor,  but  also  by  (i)  source  radiation  pattern,  (ii)  focusing/defocusing  and  (iii)  sta¬ 
tion  site  effects.  The  multiple  higher  modes  comprising  the  Lg  phase  have  radiation  pat¬ 
terns  which  might  differ  from  one  another,  but  the  gross  Lg  variation  pattern  is  expected  to 
be  moderate  (Alexander,  1985);  this  is  particularly  true  for  explosive  sources  which  can  be 
largely  represented  by  isotropic  moment  tensors.  The  station  site  effect  on  Lg  is  expected 
to  be  minimal  for  modem  broad-band  stations,  since  they  are  usually  installed  on  hard-rock 
sites,  often  in  tunnels  or  caves.  Focusing/defocusing,  on  the  other  hand,  is  more  difficult  to 
avoid  and  may  produce  apparent  Lg  Q  values  which  differ  systematically  from  values 
which  would  be  produced  by  intrinsic  anelastidty  or  by  random  scattering. 

3.  INVERSE  METHOD 

We  assume  that  we  have  a  number  of  stations,  N,  recording  Lg  from  a  single  source, 
and  that  for  the  ith  station,  a  number,  /(i),  of  Lg  spectral  estimates  are  available.  As  with 
previous  spectral  inversions  for  path-invariant  Q  (Campillo  et  al.,  1985;  Sereno  et  a/.,  1988), 
we  introduce  a  residual  square  as 

i-N  /  •/«) 

*  2  2  €2  ,  (11) 

i»l  /*1 

where  €1;  is  defined  as 


f  .  1--H, 

[^wpo/G*1^  )*;(/,)]  -  In 

S(/j)exp 

*/>  Tf 

Qa 

f  j  is  the  jth  discrete  frequency  and  A,(/p  is  the  Fourier  spectrum  for  the  windowed  time 
series  containing  Lg.  It  is  different  from  Alt,(/; )  in  that  ambient  noise  and  signals  from 
other  phases,  such  as  SR  coda,  may  be  present  in  A,(/;).  «l(  is  thus  not  identical  to 
log (ri(//))  equation  (1).  Qw  and  -q,  are  now  used  to  denote  apparent  Q  parameters  (i.e.. 
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from  nov.  on  the  prime  signs  will  be  dropped  for  simplicity).  All  of  the  remaining  symbols 
in  equation  (12)  are  the  same  as  those  defined  in  the  last  section. 

We  introduce  a  model  parameter  vector,  m: 

m  *  (M*  fc,  Qqi,  Tlt»  •••»  Qqm«  "Hn)  '  (13) 

where  the  superscript  "T"  denotes  transpose.  Note  that  equation  (13)  implies  that  for 

explosive  sources  we  assumed  a  known  3  (equation  (5)).  The  residual  defined  in  equation 

(11)  and  (12)  is  obviously  a  function  of  m,  and  the  inverse  problem  is  defined  so  as  to  find 

an  in  the  model  space  such  that 

R«2(mmin)  =  min  .  (14) 

The  inverse  problem  thus  defined  has  two  important  characteristics: 

(1)  It  is  a  highly  overdetermined  problem.  The  dimension  ofmis2xN  +  2,  N  being  of 
the  order  of  10°  to  101,  whereas  2/(i)  is  usually  in  the  order  of  103  (e.g.,  for  the  JVE 

i 

event  described  in  the  application  section  of  this  paper,  N  is  12  and  the  2/(i)  is 

i 

1822). 

(2)  It  is  a  non-linear  problem  (equation  (8),  assuming  that  A,(f)  is  reasonably  dose  to 
Aj{f  )),  and  therefore  the  residual  (equation  (11))  may  have  many  local  minima. 

The  first  characteristics  further  rationalizes  our  parameterization  of  allowing  Qu  and 
Hi  to  be  variable  among  paths.  The  second  characteristics  suggests  that  a  linearized  search 
for  will  be  successful  only  when  a  good  knowledge  on  m^,  and  thereby  a  good  start¬ 
ing  model  that  is  dose  to  ^  is  available.  This  can  be  difficult  to  obtain  when  in  addition 
to  the  variables  M0  and  fc,  Qa  and  Vi  are  allowed  to  be  variable  among  paths.  A  full  non¬ 
linear,  exhaustive  search  for  m^,  on  the  other  hand,  will  be  very  time-consuming.  We 
therefore  limit  the  non-linear  search  only  to  a  small  portion  of  m.  We  partition  m  in  the  fol¬ 
lowing  way: 

m-(nv  “V)7  - 


where 


(15) 
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*  (Mo ‘fc)T 

and 


(16) 


®p  *  (Qoi,  'll*  •••»  Qon>  ^n) 

We  now  introduce  a  quantity  A  V,  defined  as 


(17) 


A',j  •  ** 


.  wr,  U(/)G(Ai) 


(18) 


From  equation  (8)  and  assuming  once  again  that  A~*(j ;)  is  close  to  A{(f  p,  we  have 


A '«/  *  (l-/ni)ln/;-ln  (q<*  )  +  ej;-  ,  i  »1,  2  ,3  - - 1  (19) 

where  e{j  is  a  small  residual.  Generally  there  is  no  explicit  relation  between  etj  and  tijt 

/( 0  /<■) 

although  to  the  first  order  ]£e2  is  proportional  to  for  each  of  the  i  stations.  Therefore 

/* ‘  /*> 

once  mt  is  known,  approximate  values  of  the  elements  of  mp  can  be  found  by  a  linear 
regression  using  equations  (18)  and  (19).  Each  linear  regression  involves  quantities  from  a 
single,  say  the  ith,  station.  A  complete  search  for  m  can  therefore  be  conducted  using  the 
following  procedure: 


(1)  Take  a  possible  combination  of  M0  and  fc,  and  calculate  S(J)  and  A'(;  for  each  of  the 
stations  using  equations  (2),  (9)  (or  (10))  and  (18). 

(2)  For  each  of  the  stations,  linearly  estimate  Qa  and  using  A and  equation  (19). 

T  T 

(3)  Estimate  the  residual  for  m  =(m^mp)  ,  m,  being  assumed  in  step  (1)  and  being 
obtained  in  step  (2). 

(4)  Take  other  possible  M0  and  ft  values,  and  repeat  steps  (1)  through  (3)  in  a  loop,  until 
all  of  the  possible  M0  and  fc  values  are  exhaustively  searched  and  the  corresponding 
residual  values  obtained. 

(5)  Find  m^  based  on  the  comparison  of  all  residual  square  values. 

In  practice,  the  search  for  m,  can  be  accelerated  by  first  using  coarse  intervals  AM0  and 
A ft  for  M0  and  fc  defined  within  possible  ranges  and  calculate  the  residual  values.  The 
resulting  residual  square  values  (or  1-Res2  values)  are  then  displayed  (e.g.,  Section  5)  for 
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roughly  locating  the  m,  portion  of  m^.  A  refined  search  around  this  rough  location,  with 
much  smaller  intervals  AMQ  and  can  then  be  conducted  to  find  m^  to  a  realistic  preci¬ 
sion. 

Once  the  global  minimum  m^  is  obtained,  the  model  covariance  matrix,  Cm.  can  be 
obtained  using  a  linear  perturbation  suggested  by  Sereno  et  al.  (1988): 

C"  *  ~N~~N  ^rG  +  GTG  ^rG  +  XI^  '  (20) 

where  Nd  *  2/(0  is  the  total  number  of  spectral  estimates  and  Np  is  the  dimension  of  m.  G 

i 

is  the  Frechet  kernel  defined  as 

G<7*  s~,n(A(/y))  •  (21) 

omk 

The  square  roots  of  the  diagonal  elements  of  Cm  give  the  estimates  for  model  standard  errors. 
In  this  study,  the  value  of  the  damping  coefficient,  X,  was  set  to  be  0.001  times  the  minimum 
of  diagonal  elements  of  GTG.  In  the  following  sections,  I  will  present  an  application  of  the 
inverse  method. 

4.  DATA 

The  inverse  method  is  applied  to  Lg  generated  by  three  underground  nuclear  explosions 
in  eastern  Kazakhstan,  near  Semipalatinsk.  Table  1  lists  the  event  parameters.  The  largest 
event  in  Table  1  is  a  Joint  Verification  Experiment  (JVE)  event,  which  generated  Lg  signal 
with  high  signal/noise  (S/N)  ratios  at  five  broad-band  IRIS  or  CDSN  stations.  For  the  remain¬ 
ing  two  events,  only  three  stations  were  triggered  by  the  signals.  The  corresponding  great 
circle  paths  are  plotted  in  Fig.  1. 

The  windowing  and  processing  of  Lg  is  conducted  in  much  the  same  way  as  that 
described  by  Xie  and  Mitchell  (1990).  In  particular,  the  Lg  wave  train  is  isolated  by  a  20  per¬ 
cent  cosine  taper  window  with  comers  located  around  arrivals  corresponding  to  group  veloci¬ 
ties  of  3.10  km/s  and  3.60  km/s.  Fig.  2  shows  an  example  of  the  time  series  containing  Lg, 
and  the  window  used,  together  with  the  Fourier  spectra  of  the  windowed  signal  and  noise 
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prior  to  the  P  wave  arrival. 

Each  Lg  power  spectrum  is  subtracted  from  that  averaged  for  ten  noise  windows  prior  to 
P,  and  only  those  Lg  spectral  estimates  yielding  energy  S/N  ratio  higher  than  about  2  are  used 
for  further  analysis.  The  frequency  bands,  in  which  such  desirable  Lg  spectral  estimates  are 
available,  varies  between  about  0.2  to  0.35  Hz  at  their  lower  ends  and  between  about  2  to  5 
Hz  at  their  higher  ends. 

5.  RESULTS 

The  inverse  method  is  applied  to  Lg  spectral  estimates  for  each  event,  with  searches 
conducted  exhaustively  in  the  m,  subspace  as  described  earlier.  Two  successive  searches,  the 
second  one  being  refined  from  the  first,  are  usually  sufficient  to  find  with  realistic  preci¬ 
sion  (Table  1).  A  single  search  over  40  by  40  possible  source  parameter  combinations,  m^, 
takes  less  than  a  few  minutes  of  CPU  time  on  a  SUN  SPARC  1  work  station.  Fig.  3  displays 
the  result  of  one  search  in  a  30  view  of  1  -  Res2  versus  M0  and  fe .  The  value  of  P  is  set  at 
0.75,  corresponding  to  a  Poisson  medium.  Forty  three  M„  and  thirty  nine  fc  values  are 
searched,  starting  at  1  x  1022  dyne-cm  and  0.4  Hz,  with  increments  of  1  dyne-cm  and  0.2  Hz, 
respectively.  The  Af0  and  fc  values  that  are  lower  than  those  searched  lead  to  negative  Q 
values  and  are  thus  non-physical.  The  surface  in  Fig.  3  is  fairly  simple  since  the  multiple 
variables  Qa  and  n,  are  not  displayed,  yet  it  is  evident  that  convergence  to  a  local  extreme 
would  occur  if  one  uses  a  linearized  search  with  a  starting  model  located  in  the  lower  left  por¬ 
tion  (i.e.,  a  starting  model  with  a  large  fe  and  small  Mq).  When  all  of  the  unknowns  are 
taken  into  account,  the  number  of  local  extrema  is  expected  to  increase  in  a  complex  manner. 

After  numerically  examining  the  result  of  the  first  search,  a  second  search  around 
M0  *  14X1022  dyne-cm  and  fe  «  0.6  Hz  is  conducted.  The  optimized  value  of  m^  leads  to 
synthetic  Lg  spectra  that  are,  over  the  entire  frequency  band,  closest  to  the  observed  (Fig.  4). 

To  test  the  robustness  of  the  resulting  to  the  value  of  p  used,  we  also  conducted 
inversions  for  all  events  with  3  values  of  0.6  and  1.0,  respectively.  The  resulting  M0  varied 
by  less  than  10%  from  those  in  Table  1,  with  the  higher  P  value  corresponding  to  lower  M0 
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values.  The  fe  values,  when  scaled  by  ^0,  varied  by  less  than  0.01  Hz.  This  suggests  that 
the  results  of  this  inversion  are  quite  insensitive  to  varying  0  values,  and  that  an  inversion 
which  includes  0  as  an  unknown  will  unlikely  resolve  it  well. 

Comparison  of  Source  Parameters  with  Pre-existing  Information 

The  simultaneous  measurements  of  M0  and  fc  values  (Table  1)  are  the  first  ever 
obtained  using  Lg  from  underground  nuclear  explosions.  We  may  compare  these  with 
relevant  pre-existing  knowledge  to  check  for  consistency.  Assuming  that  the  M0  values 
derived  from  Lg  are  equivalent  to  those  from  other  body/surface  waves,  and  assuming  a  por¬ 
tability  of  the  M0-yield  relationship  developed  for  the  NTS  by  Denny  and  Johnson  (1991,  Fig. 
8),  the  M0  value  for  die  JVE  event  in  Table  1  corresponds  to  a  seismic  yield  of  (130  ±  10)  kt. 
This  agrees,  within  the  estimated  uncertainties,  with  the  previous  estimates  obtained  with 
seismic  and  non-seismic  methods.  For  example,  Priestley  et  al.  (1990)  obtained,  by  measuring 
time-domain  amplitudes  of  Lg  generated  by  the  same  JVE  event,  but  recorded  at  stations  dif¬ 
ferent  from  those  in  this  study,  seismic  yields  of  148  kt  and  118  kt  based  on  the  -yield 
relationship  of  Nuttli  (1986a)  and  that  of  Patton  (1988),  respectively.  The  uncertainty  of  these 
estimates  should  be  at  least  27%  (Patton,  1988).  Sykes  and  Ekstrom  (1989),  jointly  using  mb 
and  Ms,  obtained  a  yield  of  113  kt  for  the  same  JVE  event,  with  an  uncertainty  factor  of  1.28. 
Hydrodynamic  measurements  by  US  and  Soviet  scientists  for  the  same  event  yielded  values 
of  115  and  122  kt,  respectively  (cf.,  Sykes  and  Ekstrom,  1989). 

The  /c  value  for  the  JVE  events  in  Table  1,  if  scaled  by  v/0,  is  very  close  to  the  value  of 
0.7  Hz  estimated  from  Pn/Pg  waves  by  Priestley  et  al.  (1990).  A  comer  frequency  of  0.7  Hz  is 
also  consistent  with  the  moment-comer  frequency  relationship  by  Denny  and  Johnson  (1991, 
Fig.  15),  in  which  a  yield  of  120-140  kt,  for  example,  corresponds  to  a  comer  frequency  of  0.7- 
0.9  Hz. 

Since  Lg  from  the  two  smaller  events  triggered  only  3  stations,  with  lower  5/N  ratios, 
our  confidence  on  the  parameters  estimated  for  these  events  are  somewhat  lower.  The  aver¬ 
age  Qo  value  from  the  event  of  12  November  is  also  slightly  lower  than  those  from  the  other 
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two  events,  although  they  do  agree  within  the  uncertainties.  Nevertheless,  if  we  use  the 
moment-yield  relationship  of  Denny  and  Johnson  as  before,  we  can  predict,  for  the  two 
smaller  events  using  the  M0  values  in  Table  1,  yields  of  (26  ±  3)  and  (18  ±  2)  kt,  respectively. 
If  we  assume  a  minor  uncertainty  of  0.2  units  for  the  mk  values  in  Table  1,  and  use  the  mk- 
yield  relationship  by  Sykes  (Nuttli,  1986b),  the  probable  yields  will  be  between  9  and  23.7  kt, 
a  range  overlapping  our  predictions  when  the  uncertainties  are  taken  into  account. 

Finally,  we  note  that  the  fc  values  in  Table  1  decrease,  with  increasing  M0  values,  more 
slowly  than  that  predicted  by  a  scaling  relationship  of  M0  *  //,  which  was  used  previously 
for  explosions  in  Scandinavia  (Sereno  et  al.,  1988).  We  feel  that  more  work  needs  to  be  done 
on  the  Lg  moment-comer  frequency  scaling  for  both  explosions  and  earthquakes. 

Comparison  of  Q0  and  rj  Values  with  Previous  Information 

The  average  Q„  and  t|  values  in  Table  1  agree  with  one  another  within  the  uncertainties, 
and  agree  with  those  by  Given  et  al.  (1988),  Sereno  (1990)  and  Chun  et  al.  (1992),  all  obtained 
for  Lg  Q  in  roughly  the  same  area.  There  are,  however,  significant  variations  of  the  Q0  and  tj 
values  from  path  to  path  (Fig.  1).  The  Q0  values  are  above  750  to  station  OBN  and  ARU,  516 
to  HIA,  591  to  WMQ,  and  475  to  GAR.  The  uncertainties  in  Fig.  1  range  between  about  10% 
to  20%  for  Q0  values,  and  are  taken  from  the  maximum  of  the  (i)  uncertainty  estimated  using 
the  square-root  of  the  corresponding  diagonal  element  in  equation  (20),  and  (ii)  the  standard 
error  obtained  in  the  linear  regression  (equation  (19)).  Even  with  the  uncertainties  thus 
estimated,  it  seems  certain  that  the  Q0  value  for  the  path  to  station  GAR  is  lower  than  those 
for  paths  to  OBN  and  ARU.  These  correlate  with  the  regional  tectonics  since  the  former  path 
passes  through  the  Tienshan  Mountains,  whereas  the  latter  paths  lie  within  the  shield-like 
Kazakhstanian  and  East  European  plates  which  have  both  been  stable  since  Palaeozoic  time 
( Uetal .,  1982). 

Previous  studies  (Singh  and  Herrmann,  1983;  Der  et  al.,  1984;  Xie  and  Nuttli,  1988;  Rya- 
boy,  1990;  Xie  and  Mitchell,  1990)  suggested  that  the  Lg  coda  Q  generally  closely  resembles 
the  when  determined  in  the  same  area.  Xie  and  Mitchell  (1991)  obtained  a  tomographic 
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map  of  the  laterally  varying  Lg  coda  Q  in  southern  Eurasia.  Fig.  5  is  a  map  showing  the  Lg 
coda  Qq,  adapted  from  Xie  and  Mitchell  (1991).  The  paths  used  in  this  study  are  superposed 
on  Fig.  5  for  comparison.  Like  the  Q0  values  in  Fig.  1,  the  Lg  coda  Q0  values  in  Fig.  5  are 
also  subject  to  uncertainties  of  about  10%  to  15%.  Assuming  a  dose  resemblance  between  Lg 
coda  Q  and  Lg  Q,  we  can  predict  from  Fig.  5  path-variable  Q0  values  that  are,  taking  into 
account  the  uncertainties,  consistent  with  those  in  Fig.  1  for  all  paths.  For  example,  based  on 
Fig.  5,  one  can  predict  Q0  values  that  are  higher  than  about  700  to  stations  OBN  and  ARU, 
and  an  Q0  value  of  about  450  to  GAR.  One  can  also  predict  that  the  Q0  value  to  HAI  is  lower 
than  that  to  WMQ.  The  consistency  between  the  Lg  Q0  and  Lg  coda  Q0  values  suggests  that 
the  focusing/defocusing  effects  on  Lg,  such  as  those  predicted  by  Bostock  &  Kennett  (1990), 
have  not  significantly  affected  the  Lg  Q  values  estimated  in  this  study. 

6.  DISCUSSION 

The  results  of  applications  of  the  inverse  method  in  this  study  are  dependent  on  die 
validity  of  stochastic  modeling  of  Lg,  particularly  on  validities  of  the  geometrical  spreading 
factor  G(A)  (equation  (2))  and  the  simplified  Mueller-Murphy  source  model  (equation  (10)). 
The  resulting  fc,  Q0  and  q  value  obtained  in  this  study  are  consistent  with  a  priori 
knowledge,  and  this  consistency  can  be  used  to  support  these  validities.  Nevertheless,  these 
validities  deserve  some  further  discussions  for  the  darity  in  future  applications  of  the  inverse 
method. 

In  this  study  we  have  used  a  G(A)  of  the  form  of  equation  (2),  with  a  A0  of  100  km. 
These  were  empirically  proposed  by  Street  et  al.  (1975)  and  were  confirmed  by  the  synthetics 
by  Herrmann  and  Kijko  (1983),  conducted  in  a  velocity  model  for  the  central  U.S.  I  will 
hereafter  refer  to  that  model  as  model  CUS.  Strictly  speaking,  the  value  of  A0  is  dependent 
on  the  detailed  crustal  velocity  structure.  For  eastern  Kazakhstan  (EK),  results  from  both  DSS 
and  receiver  function  studies  (cf.  Priestley  et  al.,  1988)  suggest  crustal  velocity  models  which 
differ  from  model  CUS.  The  differences  include  that  (i)  the  Moho  is  dipping  in  EK,  with  an 
average  depth  that  is  a  few  km  greater  than  40  km,  the  depth  of  the  Moho  in  model  CUS;  (ii) 


the  average  velocities  in  the  upper  crust  in  EK  is  lower  than  in  CUS.  Difference  (i)  tends  to 
make  A0  in  EK  greater  than  that  in  CUS,  whereas  difference  (ii)  has  an  opposite  tendency. 
Overall,  we  expect  that  a  A0  of  100  km  used  in  this  study  should  not  have  led  to  significant 
biases  in  the  resulting  model  parameters  obtained.  However,  in  future  studies  it  will  be 
worthwhile  to  investigate  whether  or  not  a  A0  of  100  km  and  a  decay  rate  of  A~w  are  proper 
for  various  study  areas  (cf..  Bowman  and  Kennett,  1991). 

The  source  model  of  Mueller  and  Murphy  (1971)  which  was  originally  proposed  for 
dose-in  acceleration  measurements.  Sereno  et  td.  (1988)  were  the  first  to  adapt  that  model 
into  the  form  of  equation  (10)  in  this  paper.  Various  other  source  models  have  also  been  pre¬ 
viously  proposed  for  explosive  sources  (cf.,  Denny  and  Johnson,  1991).  It  is  possible  that  one 
of  the  other  source  models  could  fit  the  Lg  spectra  equally  well.  Moreover,  there  are  always 
uncertainties  in  the  value  of  as  mentioned  before,  and  in  the  values  of  p  and  v,  (equation 
(10)),  the  average  density  and  shear  wave  velodtv  in  the  crust.  In  future  studies,  the  non- 
uniqness  of  the  source  model  and  uncertainties  in  the  model  parameters  should  be  kept  in 
mind  and  should  be  reduced  using  a  priori  knowledge  when  possible. 

Another  uncertainty  in  this  study  is  in  the  estimate  of  noise  present  in  the  Lg  window. 
The  estimate  of  noise  in  this  study  is  based  on  noise  prior  to  the  P  arrivals.  The  high- 
frequency  cut-offs  in  this  study  range  between  about  2  and  5  hz.  At  higher  frequencies,  a 
more  conservative  estimate  of  noise  could  be  obtained  using  Sn  coda  (Shin  and  Herrmann, 
1987). 

7.  CONCLUSIONS 

A  non-linear  inverse  method  is  developed  for  jointly  inverting  for  Lg  source  spectra  and 
path  Q.  This  method  has  the  following  characteristics: 

(1)  It  allows  path  Q0  and  ti  values  to  be  variable  among  paths. 

(2)  It  is  an  event-based,  multiple  station  method,  suitable  for  real-time  processing  with 
modem  broad-band  stations. 


(3)  It  searches  the  source  parameters  non-linearly  and  Q  values  linearly,  requiring  no  start¬ 
ing  model.  The  results  of  the  inversion  thus  does  not  depend  on  a  priori  knowledge  of 
the  source  and  Q. 

(4)  The  search  algorithm  is  implemented  with  successive  stages  of  refinement,  making  it 
fast  enough  to  be  mini-computer  based. 

The  method  is  applied  to  three  underground  nuclear  explosions  in  eastern  Kazakhstan, 
near  Semipaladnsk.  The  resulting  seismic  moment  and  comer  frequency  for  the  largest  (JVE) 
event  are  consistent  with  available  previous  estimates  using  other  methods,  and  the  empirical 
relationships  developed  for  NTS  explosions.  Assuming  that  the  seismic  moments  obtained 
from  different  seismic  phases  are  equivalent  and  a  portability  of  various  source-scaling  rela¬ 
tionships  developed  for  the  NTS  events,  die  Lg  seismic  moments  obtained  in  this  study  for 
the  two  smaller  events  correspond  to  yields  that  are  consistent,  within  reasonable  uncertain¬ 
ties,  with  the  mb  values  by  the  ISC. 

The  path-averaged  Q0  values  and  q  values  obtained  in  this  study  are  consistent  with 
previous  estimates  (Given  et  al.,  1990;  Sereno,  1990;  Chun,  1992).  The  Q0  values  for  indivi¬ 
dual  paths  in  this  study  show  a  correlation  with  the  regional  tectonics,  and  agree  with  a  pre¬ 
vious  tomographic  map  of  Lg  coda  Q0  (Xie  and  Mitchell,  1991). 

In  future  applications  of  the  inverse  method  to  earthquake  sources,  azimuthal  variations 
of  the  source  radiation  pattern  might  be  more  significant  than  for  explosion  sources.  Accord¬ 
ingly,  interpretation  of  the  path-variable  apparent  Lg  Q  values  needs  to  be  approached  with 
caution. 

In  view  of  the  growing  high-quality  broad-band  seismic  data  base,  future  applications  of 
the  inverse  method  developed  in  this  study  will  likely  yield  source  parameters  for  increasing 
numbers  of  intraplate  earthquakes  and  explosions,  as  well  as  Lg  Q  measurements  along  vari¬ 
ous  paths.  These  should  provide  important  input  to  our  understanding  of  the  excitation  of  Lg 
by  various  types  of  sources,  at  various  depths,  and  to  our  interpretation  of  Q  in  terms  of  the 
heterogeneity,  composition  and  evolution  of  the  continental  lithosphere. 
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Table  1.  Event  and  path  parameters  t 


mb 

Location 

Seismic 

Moment  (dyne. cm) 

Comer 

Frequency 

Average 

.  On 

Average 

■n 

Sep.  14,  88* 

3H59M57.6S 

6.1 

49.81 'N,  78.80'E 

1.3  (±  0.1)  x  1023 

0.56  =  0.02  hz 

657  £  129 

0.31  r  0.02 

Nov.  12,  88 

3H30M30.8S 

5.4 

50.05'N,  78.98‘E 

2.6  (±  0.3)  x  1022 

0.82  *  0.02  hz 

557  st  94 

0.48  s:  0.01 

Nov.  23,  88 

3H57M06.8S 

5.4 

49.77*N,  78.06*E 

1.8  (±  0.2)  x  1022 

0.70  -  0.02  hz 

621  ±  127 

0.49  ±  0.02 

t  The  origin  times,  locations  and  magnitudes  are  from  the  ISC  bulletin;  the  seismic  moments,  comer  frequencies  and 
path-averaged  Qq  and  11  values  are  obtained  in  this  study  with  (3  in  equation  (5)  set  to  be  0.75  (Poisson  media). 

*  The  Joint  Verification  Experiment  (JVE)  event. 
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Figure  Captions 

Fig.  1.  Great  circle  paths  from  the  three  nuclear  explosions  to  the  IRIS/CDSN  stations. 
Numbers  near  the  paths  are  the  Qq,  and  q,  values  obtained  for  the  path,  with  the  estimated 
uncertainties. 

Fig.  2.  (Top)  Time  series  containing  Lg  from  the  JVE  event  observed  at  IRIS  station  ARU.  The 
trace  is  instrument  uncorrected,  and  starts  at  group  velocity  of  4.1  km/s,  corresponding  to  a 
lapse  time  of  374.0  sec.  The  smooth  curve  represents  a  20  percent  cosine  window  used  in  the 
analysis.  (Bottom)  The  instrument-  uncorrected  Fourier  amplitude  spectrum  of  the  signal  con¬ 
taining  Lg,  and  the  Fourier  spectra  from  several  noise  windows  prior  to  P. 

Fi.  3.  1  -  Res2  versus  M0  and  fc  resulting  from  a  search  over  4J  M0  values  (between  8  to  50 
x  1022  dyne-cm)  and  39  fc  values  (between  0.4  to  8.0  Ha),  with  respective  intervals  of  1  xio22 
dyne-cm  and  0.2  Hz.  The  £  value  (equation  (5))  is  set  to  be  0.75.  Note  that  the  Qa  and  q, 
values  are  not  displayed,  making  the  surface  considerably  simpler,  with  fewer  local  extrema. 

Fig.  4.  Synthetic  Lg  spectra  for  all  of  the  five  stations  recording  the  JVE  event,  versus  the 
observed.  The  spectra  are  normalized  to  source  by  multiplying  by  4-rrpo,  (ie.,  argu¬ 

ments  of  the  first  and  second  logarithms  in  equation  (12)  are  the  observed  and  synthetic  spec¬ 
tra  displayed  here,  respectively).  The  units  used  are  (Hz)  for  frequency  and  (1022  dyne-cm) 
for  amplitude.  The  synthetic  spectra  are  calculated  using  the  optimal  source  model  parame¬ 
ters  in  Table  1,  and  optimal  Q^,  q,  parameters  indicated  above  each  panel. 

Fig.  5.  Tomographic  map  of  laterally  varying  Lg  coda  Qb  in  the  area  under  study,  adapted 
from  Xie  and  Mitchell  (1991).  The  paths  used  in  this  study  are  superposed  for  comparison. 
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